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Approaching a dangerous bifurcation, from which a dynamical system such as the Earth's 
climate will jump (tip) to a different state, the current stable state lies within a shrinking 
basin of attraction. Persistence of the state becomes increasingly precarious in the presence 
of noisy disturbances. We argue that one needs to extract information about the nonlin- 
ear features (a "softening") of the underlying potential from the time series to judge the 
probability and timing of tipping. This analysis is the logical next step if one has detected 
a decrease of the linear decay rate. If there is no discemable trend in the linear analysis, 
nonlinear softening is even more important in showing the proximity to tipping. 

After extensive normal form caUbration studies, we check two geological time series 
from paleo-climate tipping events for softening of the underlying well. For the ending of 
the last ice age, where we find no convincing linear precursor, we identify a statistically 
significant nonlinear softening towards increasing temperature. 

Keywords: time series analysis, bifurcation prediction, climate tipping 

1. Introduction 

The on-going work by the United Nations, following the major report of the IPCC [1] in 
2007 and the Climate Change Conferences in Copenhagen and Cancun (Mexico), centres 
on the prediction of future climate change, a key feature of which would be any suggestion 
of a sudden and (possibly) irreversible abrupt change called a tipping point (Lenton et al. 
[2], Scheffer [3]). Many tipping points, such as the switching on and off of ice-ages, are 
well documented in paleo-climate studies over millions of years of the Earth's history. 

Predicting such tippings in advance using time-series data derived from preceding be- 
haviour is now seen as a major challenge, impinging, for example on the possible use of 
geo-engineering (Launder and Thompson [4]). Techniques introduced by Held and Kleinen 
[5] and Livina and Lenton [6] draw on the assumption that the tipping events are governed 
by a bifurcation in an underlying nonlinear dissipative dynamical system. Specifically, re- 
searchers search for a slowing down of intrinsic transient responses within the data, which 
is predicted to occur before most bifurcational instabilities. This is done by determining 
a so-called propagator, estimated from the correlation between successive elements of a 
window sUding along the time series (this estimate is called jcacf in this paper). This prop- 
agator, such as the AR(1) mapping coefficient, is a measure of the linear decay rate and 
should increase to unity at a tipping instability (corresponding to a decrease of the linear 
decay rate to zero). 

Prediction techniques can be tested on climatic computer models, but more challeng- 
ing is to try to predict real ancient climate tippings, using their preceding geological data 
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provided by ice cores, sediments, etc. Using this preceding data alone, the aim would be to 
see to what extent the actual tipping could have been predicted in advance. 

One such study by Livina and Lenton [6] looks at the end of the most recent ice age 
and the associated Younger Dryas event, about 11, 500 years ago, when the Arctic warmed 
by 7°C in 50 years. This pioneering study used a time series derived from Greenland ice- 
core paleo-temperature data. A second such study (one of eight made by Dakos et al. 
[7]), using data from tropical Pacific sediments, gives a good prediction for the end of 
'greenhouse' Earth about 34 million years ago when the climate tipped from a tropical state 
into an icehouse state. [Lenton et al. this volume] gives a complete overview of current 
techniques for extracting early-warning signals from time series and compares them on 
realistic models and a range of paleo-cUmate time series. 

A remark that is made during the analysis of [Lenton et al. this volume] is that the 
absolute values of the extracted quantities have no direct bearing on the tipping probability 
over time unless one knows also the size of the basin of attraction. This basin is determined 
in the simplest cases by the dominant nonlinear term in the underlying equations of motion, 
which is the subject of this paper. 

2. Concepts from Nonlinear Dynamics 

The core techniques for analysis of time series of cUmate data aim to extract the linear de- 
cay rate toward (a possibly drifting) noise-perturbed equilibrium [5, 6]. These techniques 
can be directed at the Earth's climate as a whole, or at the relevant cUmate sub-systems 
described by Lenton et al. [2] as tipping elements. Our aim here is to augment this linear 
information with information about nonlinear features of the underlying dynamics, by ex- 
amining to what extent we can identify nonUnear softening and include it into the hst of 
early-warning signs for climate tipping. 

(a) Shrinking basins around dangerous bifurcations 

As we approach a dangerous bifurcation, from which a nonlinear dissipative dynamical 
system will experience a finite jump to a remote alternative state, the current attractor is 
located within a shrinking basin of attraction. Correspondingly, the maintenance of the 
state will become increasingly precarious in the presence of noisy disturbances. 

Thinking in terms of an underlying potential energy function (the existence of which is 
theoretically well-defined for a saddle-node fold and some other bifurcations) the parabolic 
shape of its graph will become increasingly perturbed by softening nonlinear features. 

The relevant bifurcations are the so-called codimension-one events (Thompson et al. 
[9], Thompson and Stewart [10]), namely those that can be typically encountered un- 
der a gradual change of a single control parameter. The complete list of the dangerous 
codimension-one bifurcations is given in Table 1, where we give a brief description of 
the nature of the shrinking basin. More details can be found in Thompson and Sieber 
[11]. Focusing first on the local bifurcations (headings (a) and (b) of Table 1) the shrinking 
boundaries are illustrated schematically in Figure 1, with the control parameter plotted hor- 
izontally and the response plotted vertically. In the first picture. Figure 1(a), we illustrate 
both types of the saddle-node fold bifurcation: the static fold involving a path of equi- 
hbrium fixed-points, and the cyclic fold involving a trace of periodic orbits. For these the 
basin shrinkage is one-sided. The basin is bounded by the unstable side of the parabolically 
folding path. 
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(a) Local Saddle-Node Bifurcations 

Static Fold (saddle-node of fixed point) one-sided basin shrinkage 

Cyclic Fold (saddle-node of cycle) one-sided basin shrinkage 

(b) Local Subcritical Bifurcations 

Subcritical Hopf complete basin shrinkage 

Subcritical Neimark-Sacker (secondary Hopf) complete basin shrinkage 

Subcritical Flip (period-doubling) complete basin shrinkage 

(c) Global Bifurcations 

Saddle Connection (homoclinic connection) outside shrinkage around cycle 

Regular-Saddle Catastrophe (boundary crisis) fractal basin shrinkage 

Chaotic-Saddle Catastrophe (boundary crisis) fractal basin shrinkage 

Table 1: Dangerous Bifurcations and the behaviour of the basin of attraction 
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Figure 1 : Basin boundary transformations at dangerous bifurcations, (a) the fold bifurcations and (b) the subcrit- 
ical bifurcations. Solid curves denote stable paths while broken curves denote unstable paths. The basin is shown 
in grey. 



Next, in Figure 1(b), we illustrate the local subcritical bifurcations of which the Hopf, 
flip and Neimark are codimension-one (we remember that the more familiar static pitchfork 
bifurcations are codimension-two, being structurally unstable against a symmetry-breaking 
perturbation). Here we see an unstable solution, which controls the basin boundary, shrink- 
ing parabolically around the stable solution, from which the (noise-free) system jumps out 
of our field of view as the control parameter, ^, reaches its critical value of Hc- 

The static fold and the Hopf bifurcation have a theoretically well-defined potential 
energy surface that neatly summarizes the shrinking basin as illustrated in Figure 2. Note 
that Figure 2(a) will be discussed more fully in section 2(b). 

{b) A closer look at the fold 

The fold is the simplest and most common way in which an equilibrium of a nonlinear 
dissipative dynamical system can lose its stability. Having only a single active degree of 
freedom (x) and being observable under the variation of a single control parameter {jx) 
it can be illustrated on a graph of x against ;U as a smooth curve that simply reaches an 
extreme value (a maximum, say) of the control jU, as shown in Figure 1(a). So although it is 
traditionally called a saddle-node bifurcation, it does not exhibit any obvious 'bifurcation' 
of a path, but rather just a smooth folding of an existing path. In the case of a static fold (to 
which we shall largely restrict our attention) the path is a trace of equilibrium fixed points, 
while for a cyclic fold the path would be a trace of periodic solutions. 

Assuming that we are close to a fold, the intrinsic damping of the system will have be- 
come super-critical so the system will have the non-osciUatory response of an over-damped 
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Figure 2: Total potential energy transformations in (a) the saddle-node fold and (b) the Hopf bifurcation (via 
averaging). Black balls denote stable equUibrium states while white balls denote unstable states. 

particle sliding (in thick oil, as we might imagine) on a notional potential energy surface 
U{x,n). This is illustrated for the fold in Figure 2(a), where we show the potential en- 
ergy surface erected over the (x, ji) base plane. Here the solid part of the equiUbrium path, 
corresponding to energy minima, is a stable node, while the dashed part, corresponding to 
energy maxima (more generally a geometrical saddle in higher dimensions) is an dynami- 
cally unstable saddle. 

We see that as the control is slowly increased towards its critical value, Hc, the stable 
equilibrium solution gets closer and closer to the hill-top potential barrier, and so gets 
progressively more precarious in the presence of noisy disturbances. With ji increasing at 
an infinitesimal rate, noise will ensure escape over the hill- top before n reaches Hc- If, 
however, ji increases rapidly, this early escape may be delayed or supressed as we have 
examined quantitatively elsewhere (Thompson and Sieber [12]). 

To give a greater perspective to our view of the fold it is useful to introduce another con- 
trol parameter jj.2 such that our parameter space is now represented jointly by the (p.\ ,112)- 
plane. We can now erect an equilibrium surface, of height x, over the two-dimensional con- 
trol plane as illustrated in Figure 3. Here we see two fold lines whose projections divide 
the control plane into regimes exhibiting either one or three co-existing solutions. These 
fold lines coalesce and vanish at the cusp. This cusp is a codimension-two phenomenon, 
which is typically only observed when we have independent control of both parameters 
(jUi and 112)- This is made transparent by the fact that a typical path in control space cannot 
be expected to pass precisely under the cusp point; by contrast, if a parameter path passes 
through a fold line all near-by paths also cross this fold line, which is what 'the fold has 
codimension-one' means. For an appropriate choice of path in the (jii , /i2)-plaiie through 
the cusp one sees a pitchfork bifurcation (Thompson and Stewart [10]). 

We will consider the case in which the main equilibrium sheet (in grey) is stable while 
the narrow inverted sheet (white) is unstable. 

A cUmate system will have a large number of possible control parameters, but in trying 
to predict a tipping point we will be studying a particular evolution which corresponds to 
a particular route in control space; this is precisely why we are only likely to encounter a 
fold, and not a cusp in Figure 3. To explore possible scenarios that might underlie a tipping 
study, it is now useful to exanoine some possible routes in our two-parameter model of 
Figure 3. 
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Figure 3: The cusp catastrophe and its associated pattern of folds. The notional axes show how such a cusp might 
arise in a climate tipping scenario. 

Consider, first, the simple scenario in which the system parameters take path BB lead- 
ing transversally through a fold line. This is the classical encounter with a fold that is 
implicitly envisaged in earlier work. In the presence of noise a premature tip from N oc- 
curs with a certain probability depending on the noise level and the speed with which the 
path is traversed, as analysed by Thompson and Sieber [11]. If the noise level were partic- 
ularly high, a jump from N could easily be perceived as a purely noise-induced jump with 
no bifurcation, even though it is actually caused by the proximity to the fold: this could 
hold even if there were no movement in the control space at all, the system having rested 
at N throughout ([Ashwin et al. this volume] call this N-tipping). Another scenario, not 
specifically illustrated in the figure, could arise if the route in control space approached 
the fold Une, but then turned back away from it. Analysis would then show a temporary 
decrease of the strength of the attraction (with the AR(1) coefficient moving towards unity) 
which might be discounted as a false alarm, even though the system did approach a fold 
and a noise-induced escape was temporarily probable. 

Finally, in the scenario of path AA, one would observe a temporary increase of the 
AR(1) coefficient, even though no bifurcation is apparent but rather a gradual shift of the 
steady state. This shows that a rapid transition can be related to a nearby fold that was 
just barely missed. Similarly, the scenario called R-tipping by [Ashwin et al. this volume] 
occurs if one changes the other control parameter in Figure 3 sufficiently rapidly from the 
path AA to N (avoiding the folds). In this scenario the system jumps to the other stable 
equihbrium despite having never crossed a bifurcation curve. 

3. Time series analysis for leading order terms — qualitative overview 

The extraction of the underlying dynamics from time series is divided into subtasks of in- 
creasing difficulty. Assume that we have a time series coming from a process that is either 
stationary or has a slowly drifting underlying parameter as suggested by the paths in 
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Figure 3. In order to predict (or identify) the dynamics one attempts to extract the leading 
orders of the dominant components of the underlying equations of motion. We illustrate 
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(b) Nonlinear times series, 
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Figure 4: Comparison of time series generated by linear (left column, (a-1) to (a-4)) and nonlinear process (right 
column, (b-1) to (b-4)). Both time series are stationary and have identical linear decay rates. As measures for 
nonlinearity we used skewness and the coefficient Cemp (as defined below in equation (3.5)). Parameters: /i(0) = 2, 
N = 1 223, A? = 0.1, fT = 1, window length for linear and nonlinear analysis, w = N/2, e = 0.01. 



the steps using a time series obtained from a process which has the saddle-node normal 
form with a slowly drifting normal form parameter jj. as its deterministic part (correspond- 
ing to Figure 2(a)) and is influenced by additive Gaussian noise. We permit the parameter 
H to drift slowly with speed e: 



dx -■ 
dM 



[;U — jc^] df + adWt, where 



edt. 



(3.1) 



Here x is the response, ji is the control parameter, and t is the time. The Gaussian perturba- 
tion dWt has zero mean and unit variance such that a controls the noise induced variance 
added to the dynamics of the deterministic saddle-node form x — ji —x^, which has the 
stable equilibrium jCgq — V/I- The second equation prescribes the linear sweep of ji from 
its starting value flo>0 towards zero at rate e <C 1. For fixed > the unperturbed system 
(C7 = 0) has a stable equilibrium atx— x^q — \/lJ. and an unstable equilibrium at x = — ^/Ji. 

In Section 4 we will use the saddle-node normal form (3.1) to test the ability of a range 
of estimators to distinguish a time series generated by a process with a quadratic nonlin- 
earity from a linear time series. Two example time series and their analysis are shown in 
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Figure 4. Both time series have identical hnear properties but differ in the dominant non- 
linear term of the deterministic part of the dynamics.The softening of the full (estimated) 
nonlinear potential well (corresponding to the illustration in Figure 2(a)) is displayed in 
Figure 5. 

(a) Equilibrium — zero order 

The location of the slowly drifting equilibrium and its trend is estimated first. [Lenton 
et al. this volume] call this step detrending and propose either piecewise linear fitting or 
filtering with a Gaussian kernel. The result of filtering with a Gaussian kemel is shown 
as a thick grey curve in the row 1 of Figure 4. Every method for detrending requires a 
bandwidth parameter (for example the width of the Gaussian kernel). For Figure 4 we 
chose the bandwidth which the kemel density estimation routine (Matlab's kde by Botev 
et al. [14J) offered. The thin black curve in Figure 4(b-l) shows the true value of the quasi- 
equilibrium (for e = 0). It gives a visual estimate of how the reliability of the zero-order 
estimates depends on the quality of the data. The oscillatory deviations of the extracted 
equiUbrium (thick light grey curve) from its true value (thin black curve) shows that the 
automatically chosen bandwidth was sUghtly smaller than the theoretical optimum. 

(b) Linear decay rate — first order 

Assuming that the deterministic part of the underlying process has a stable equilibrium Xgq 
(which is possibly slowly drifting), and that the zero-order estimate is an approximation 
for Xeq, the next step is to estimate the linear decay rate K toward Xgq and the trend of K 
over time. Generally, if the underlying deterministic process is high-dimensional then these 
estimates capture the dominant (that is, smallest) decay rate. They are typically applied to 
the detrended time series, that is, x — x — x^, which is shown in the second row of Figure 4. 

• ACF The A:-step (usually one-step) autocorrelation function (ACF) a is fitted directly 
to the detrended time series: Xk+i = CCXj^. This was introduced by Held and Kleinen 
[5] as degenerate finger-printing for the analysis of climate time series. The ACF a 
is related to K via 

a = exp(-KACFAr). 
The estimate for ktacf is shown in Figure 4(a-3) and (b-3). 

• DFA Detrended fluctuation analysis has been introduced by Livina and Lenton [6] 
for cUmate time series, see also [Lenton et al. this volume] for a short description of 
how one extracts K from this analysis. 

• Quasi-stationary density If one assumes that the deterministic dynamics of the 
detrended time series x is essentially that of an overdamped particle moving in a 
potential well U (x), that is, 

djc=-a;,f/(jc)d/-h(7dW,, (3.2) 

then the Unear decay rate equals the second derivative dxxU (0) of U at the bottom of 
the well, x = 0. Moreover, the stationary density p{x) is related to the shape of the 
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potential well U via the stationary Fokker-Planck equation 

dx [—a~^dxU {x) p{x)^ , which gives, integrated once, 

-a-^dxU{x)p{x) + (y-\ (3.3) 

where c is the constant of integration. This constant satisfies c = limx^±ac,diU {x)p{x), 
and can be interpreted as the flow rate (c < indicates flow toward —0°). Remember- 
ing that the drift speed e of the system parameter ^ is small, the true time-dependent 
probability density p stiU satisfies (3.3) approximately at each /i: 

^dxp{x,n) = -a-^dxU{x,n) p{x,n) + a-^cin) + 0{e). (3.4) 

In (3.4) we included the dependence of all quantities on ^ expressly and remember 
that fi{t) is of order e. Assuming quasi- stationarity, we neglect the order-e terms 
and determine the coefficients, —a^^dxU and a^^c, by fitting the empirical density 
Pem:p{x) from a time window [t — w/2,t + w/2] to the Equation (3.4). Specifically, 
Cemp = o~^c is a scalar and dxU (0) is equal to after detrending such that one has 
to fit two coefficients. Kg and Cemp using (3.4) if one truncates dxU after first order: 

^^^Pemp(i) = -K'C/.^Pemp(.^) +Cemp. (3.5) 

In problems where no escape is possible one can drop the term c in (3.3) (and, thus, 
in (3.4) and (3.5)), and simplify (3.3) to 

U = -^\ogp. (3.6) 

Livina et al. [15] used relation (3.6) (fitting of U with higher-order polynomials) to 
detect potential wells in time series that included rapid transitions. Fitting U with a 
parabola (where U{Q) is irrelevant and dxU{0) = 0) is equivalent to using the vari- 
ance (Tgnip of Pemp{x) if the density pemp is Gaussian. Ditlevsen and Johnsen [16J 
state that monitoring the empirical variance (Jg^ip helps to avoid erroneous detection 
of false trends. 

[Lenton et al. this volume] compare these three estimates, K"acf, the DFA-based estimate 
for K, and the variance Cg^p (which is inversely proportional to Ku), in detail using time 
series arising in climate models and in palaeocUmate records. Methods based on properties 
of the spectrum of the time series were proposed and investigated by Kleinen et al. [17], 
Biggs et al. [18] but are not discussed here. Figure 4 shows in row 3 how the estimates of Ku 
from the quasi- stationary density compare to the ACF estimates ktacf- For both time series 
the estimates are quantitatively close to the theoretical value of K (which is indicated by 
the thin black fine in Figure 4(a-3,b-3)). We also observe that in Figure 4(a-3,b-3) the local 
trends of ktacf from the autocorrelation estimate and of % based on the quasi- stationary 
density are strongly correlated, so it is unclear if monitoring both quantities really prevents 
false alarms. 

There is one notable difference between the AR(1) and the DFA estimate on one hand 
and the distribution-based estimate fQ/ on the other hand: the estimate for Ku based on the 
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quasi- stationary density only looks at the distribution of the values in the window of 
interest but does not care about the order in which they appear. This is in contrast to, for 
example, the AR(1) estimate which determines the correlation between each value and 
its predecessor. 

Estimate of the input noise amplitude 

If the AR analysis shows the presence of a single positive dominant coefficient (which gives 
evidence that the underlying dynamics has really one distinct direction in which the decay 
is slowest) then one can combine both estimates of the decay rate to obtain an estimate of 
the amplitude o of the additive noise. Equations (3.3) and (3.5) imply that fQ/ is an estimate 
for a~^K where k and are the true linear decay rate and input noise amplitude. Thus, 
if we replace the unknown K by the estimate ktacf we can use the relation between Ku and 
the true K to estimate a^: 

<7emp = KaCf/Ku (3.7) 

where Kacf is the estimate of the linear decay rate obtained from the autocorrelation and 
fff/is the estimate obtained from the quasi-stationary density /?emp- This is an alternative to 
using the residual of the least-squares fit that one obtains when estimating Kacf- We found 
that the residual systematically underestimates the noise level in the normal form examples 
shown in Fig. 4. 

(c) Nonlinearity — basin boundary 

One problem left open by the linear analysis is that the quantitative value of the estimated 
linear decay rate k and the estimated noise-level do not give any indication for the proba- 
biUty of tipping. What the Unear methods estimate is the a and the C7 of a discrete process 

Xk+i = axk + arij, (3.8) 

(where the rji^ are assumed to be independent and normally distributed random numbers). 
The discrete process (3.8) is the linearisation of the time-A? map of the continuous process 
(3.2). An estimate for a less than but close to unity does not necessarily indicate that we 
are close to a stability boundary for the nonlinear problem (3.2). Rather, it implies that 
the time step At between successive measurements has been small compared to the mean 
decay time to half, log2/K', of the process: 

a{tk) = exp {-K{tk)At) = 1 - K{tk)At + O {{Atf) . 

The trend of the estimated a is an indicator for incipient tipping because (if extrapolated in 
any way) it gives an estimate for the time to tipping if one ignores the possibility of early 
escape. However, the trend of a is less certain even in the artificial time series of Figure 4. 
Similarly, the estimated noise amplitude a has to be compared to the coefficient in front 
of the leading nonlinear term of the right-hand side (which equals —1 in (3.1)). The ratio 
between o and nonlinear term enters the estimates for the probabihty of tipping if there is 
no discernible upward trend in a (N-tipping), or for the probability of early escape if the 
trend of a points upwards (see Thompson and Sieber [12]). So without at least an order- 
of-magnitude estimate of the leading nonhnear term the hnear methods only provide an 
estimate for the tipping time if there is a clearly discernible trend in a, and even then this 
estimate discounts early escape. 
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Figure 4 illustrates this problem. Apart from a slow trend both time series in Figure 4 
have identical properties at the linear level. However, while for the time series in Fig- 
ure 4(b-l) the median time to tipping is f = 100 according to Thompson and Sieber [12] 
(indeed, it escapes shortly after t = 120), for the time series in Figure 4(a-l) the extrapola- 
tion (correctly) predicts tipping (or rather, linear instability) at f « 320. Thus, we observe 
that the linear properties of the time series alone are not sufficient to provide estimates 
about tipping, for example, the median time to tipping. The difference between the time 
series in Figure 4(a-l) and Figure 4(b-l) is in the dominant nonlinear term of the deter- 
ministic part. That this term plays a role is not surprising because tipping is a feature of 
nonlinear systems. Extracting its presence from a time series is a necessary step that has 
to follow the linear analysis. The approaches presented in the following all generalise 




Figure 5: Nonlinear potential energy surfaces extracted from noisy time series using sliding windows (w = N/2) 
for the saddle node fold (f/cmp = ^f^cmp log(Pcmp) /2 where CTcmp was estimated using (3.7)). At either end of the 
surface we have plotted the known function U {x,t) = — /i(f)jc (a = 1, e = 0.01 and /i(0) = 2). Red denotes 
a positive deviation, blue a negative deviation from the parabola given by linear theory. 

the estimate for Kjj, based on the quasi-stationary density, by looking at the full nonlinear 
potential well f/emp(-^) associated to the quasi-stationary density Pemp(-^) via the Fokker- 
Planck equation (3.4). Figure 5 illustrates how t/emp = — log(/?emp) /2 looks like for the 
time series shown in Figure 4 column (b). The realisation of the noisy, evolving time series 
is shown in the base plane together with the equilibrium path ji = x^, and the estimated 
(empirical) potential, t/emp is shown as the coloured surface. The estimate for the potential 
well is taken in a sliding window, which explains why the surface can only be determined 
in the central region of the time series with half the length of the window unrepresented at 
either end. The numerically estimated equilibrium trend Xeq is displayed as a black curve at 
the bottom of the valley of t/gmp- For comparison with this t/emp we display at either end the 
real (nonlinear) potential U{x,t) corresponding to (3.1), which at the later time is already 
showing the hill-top to the left. The colour in Figure 5 shows the deviation of the empiri- 
cal potential from the linear-theory parabola, which is not itself displayed. The softening, 
highlighted by the blue colouration, is seen on the left hand side which is approaching the 
hill-top. Transparency of the colour is used to show the number of data points that support 
a particular part of the surface t/emp (■^) (as given by the empirical quasi-stationary density 

Pemp{x)). 
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Several methods have been proposed to detect nonlinearity in the potential well in the 
form of this type of softening (or tilting). 

• Skewness Guttal and Jayaprakash [19, 20J studied how the closeness of a bifurcation 
influences the skewness 7 of the empirical quasi-stationary density from time series. 
As this approach is also based on the empirical quasi- stationary density />emp. it 
generahses the estimate Ky to non-parabolic features of the well U and non-Gaussian 
features of the stationary density. Guttal and Jayaprakash [19, 20] proposed to look 
at the trend of the skewness y over time to detect incipient bifurcations. Figure 4(b- 
4) shows how the skewness, taken over a single window, changes. It is noticeably 
negative (compare to the empirical skewness from the linear time series in Fig 4(a- 
4)) but no trend is discernable. 

• Quasi-stationary density Livina et al. [15] generalised the potential well analysis 
for the quasi-stationary density beyond the estimate Ku for the linear decay rate, 
fitting the potential well to higher-order polynomials of even degree using rela- 
tion (3.6). This analysis was not used to attempt prediction of transition probabilities 
from time series but it was applied to time series that already included all transitions 
to count the number of wells of U depending on time t (which corresponds to the 
number of modes (peaks) of the empirical quasi-stationary density /?cmp)- See also 
[Beaulieu et al. this volume] for methods that detect change points in time series. 

• Drift ratio If one assumes that the underlying system parameter is approaching a 
fold more or less linearly (jL[(f) = e > 0) and one is sufficiently close to the fold 
already then the ratio between the drift A.v'cq of the equilibrium estimate and the 
increase of the estimated hnear decay rate Ajc estimates the quadratic term in the 
normal form [12]. 

4. Quantitative estimates of the nonlinear coefficients 

The colour shading of Figure 5 suggests that the nonlinear term in the underlying deter- 
ministic system shows up in the quasi- stationary density Pemp- This raises two questions. 
First, can the level of deviation from linear theory be distinguished with confidence from 
random chance? That is, can one state, when looking at a time series, that a significant 
quadratic term is present? Second, is it possible to quantify the size of the nonlinear term 
with reasonable level of certainty? 

Finding a significant nonlinearity is far easier than estimating its precise value because 
biases introduced during the analysis of the equiUbrium and decay rate (zero-order and 
first-order analysis) may distort the value of the nonlinear estimate but may still keep it 
significantly different from zero. We will demonstrate this in Section 5 for two paleo- 
chmate time series. 

(a) Properties of the stationary density of the saddle-node normal form 

Figure 6 shows the dependence of the mean, the variance varp, and the skewness 7 of 
the empirical stationary density /?emp on the system parameter /i for (3.1) as computed us- 
ing the stationary Fokker-Planck equation (3.3). Note that this corresponds to the stationary 
process where the random realisation is re-injected at -|-oo whenever it has escaped to —00. 
We observe that the dependence of 7 is strongly nonlinear for the saddle-node normal form 
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Relation between skewness and system parameter fi 





-2 

1 2 /i 3 4 

Figure 6: Dependence of mean, variance varp and skewness / of tlie empirical stationary density for the saddle- 
node normal form with additive Gaussian noise (3.1) as computed directly from the Fokker-Planck equation (3.3) 
with — (9vf/(-i;) = fl—x^ and fj = 1. The escape rate c is defined by the constant of integration in Equation (3.3). 

with noise = e = 0, CJ = 1 in Equation (3.1)). This means that trends of the skewness 
are Ukely to be difficuh to ascertain when approaching a fold. However, the presence of 
skewness is an indicator for the nonlinearity uniformly for all parameters /i, and its sign 
indicates the direction of escape. 

Note also how the nonlinearity affects the empirical mean and the variance var p, which 
is no longer inversely proportional to the linear decay rate K. This shift of the mean and 
the additional variance is in part due to the fraction of trajectories that escape, in part due 
to the non-parabolic shape of the well at some distance from its local minimum. 

ib) Practical estimates of the nonlinear term from time series 

We can estimate the nonlinearity directly as the deviation of the empirical potential 
well U from a parabola as suggested by figure 5 (we are only interested in d^U). The 
simplest approach (apart from looking at skewness 7) is to fit Ku and Cemp to the empirical 
density psmp using Equation (3.5). If dxU{x) is indeed linear then Cemp would be zero such 
that Cemp is a signed scalar measure for the deviation of dxU{x) from linearity, serving as 
a proxy for the present nonlinearity. A second alternative is to expand —d^U (i) to second 
order, incorporating a quadratic nonlinearity directly into d^U with an unknown coefficient 
A'2. This leads to 

—a^^dJJ{x) = —Kux + N2x^, such that we fit 

^d:,Pemp{x) = [-KuX+N2X^] Pemp(^) +Cemp,2- (4.1) 

Both quantities, Cemp and N2, measure the deviation of ~a^^djU{x) from its linear ap- 
proximation K[jx. The main difference between them is the weighting (x^ for N2, uniform 
for Cemp)- In summary, the step-by-step procedure to extract an estimate of the quadratic 
term (or a proxy for it) from a time series Xk is the following. 
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1. Detrend the time series xj. (see Section 3(a)). We call the detrended time series x^.. 
The following steps are performed for each suitable time t to obtain fQ/(r), Cemp(f)' 
N2{t) or skewness 7(f): 

2. Choose a length w of sUding windows and obtain the empirical stationary density 
/?emp(-^) for the sliding window centered at f . (We used kdeld to obtain /?emp(-^)- For 
non-stationary time series we used the sliding window length w — N/2.) 

3. Use relation (3.5) to find Ku{t) and Cemp(f ) from pemp- Use relation (4.1) to findN2{t) 
(and another estimate for Ku) from pemp- Compute the empirical skewness 7(f) from 

Pemp- 

(c) Comparison of estimates using saddle-node normal form 



X 10 



2 

_ 
-2^ 

0.5 

-0.5 . 
-1 
0.2" 

- 

-0.2- 

-0.4 _ 
2000 




mean 

■75%. 




Cg^p for linear surrogate 

Cpn-n for saddle-node normal form 



N2 for linear surrogate 
- for saddle-node normal form 



skewness y for linear surrogate 
- skewness y for saddle-node normal form 



length of time series before escape 



2V¥ 

Kacf (25%, 50%, 75%) 
Ku (25%, 50%, 75%) 



^1 1.5 2 fj, 2.5 3 

Figure 7: Quardles of lengths (d), and estimates c^^p (a), N2 (b), skewness 7 (c), iCacf, and Ku (botli (e)) for 
time series generated by the saddle-node normal form (3.1) with noise (e = 0, = 1) and for linear time series 
generated by (4.2). Without nonlinearity transition would occur at fl =0. 



Figure 7 shows how the estimates for nonlinearity behave for the saddle-node normal 
form with noise. Equation (3.1) with e = and C7 = 1. Each panel shows the quartiles of 
the distribution of the estimate for 100 realisatons of time series generated by (3.1). Each 
realisation was run until it reached 2000 data points or x = —1 (indicating escape). Panel 
(a) shows the estimate Cemp obtained by linear-least-squares fitting of the approximate sta- 
tionary Fokker-Planck equation (3.5) to the empirical stationary density pemp, panel (b) 
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shows the quadratic coefficient ^2 obtained from (4.1), and panel (c) shows the skewness 
7 of Pemp- For comparison, we generate also 100 linear time series using 

d;c = -l^xAt + m, (4.2) 

and then fit Cemp. N2 and 7 for these time series, too. Figure 7 compares the quartiles of 
the distributions for time series generated by the saddle-node normal form (3.1) and the 
Unear equation (4.2). If there is small overlap in the distribution then the quantity is a good 
indicator for the presence of a quadratic nonlinear term. We observe that this is the case 
for all quanitites as long as the time series has a length of 2000 (see Figure 7(d) for the 
distribution of time series lengths). Naturally, the uncertainty, and, hence, the overlap, in- 
creases when the time series is shorter due to escape from the potential well to —0° . This is 
the case for smaller ^ in Figure 7 (see Thompson and Sieber [12] for a quantitative study 
of escape probabilities for small jU). We also notice that all quantities have a systematic 
deviation from the theoretical value of the quantity they supposedly estimate: N2 should 
be equal to —1 according to (3.1), the real escape rate c has a much smaller modulus than 
Cemp (compare Figure 6(a)), and the skewness 7 has theoretically a much larger modulus 
(compare Figure 6(b)). The deviation for A^2 is relatively small, and is hkely caused by 
the kernel density estimate for Pemp- The large modulus of Cemp shows that the linear term 
— KuxK a poor fit for —dxU{x), which makes Cemp a measure of how much the odd symme- 
try is broken by —dxU (x), rather than an estimate for the escape rate. We note that Cemp,2 
as fitted in (4.1) has a reaUstic modulus (close to zero, not shown in Figure 7). The bias of 
the skewness y is mainly due to our restriction to time series that stay inside the potential 
well. However, the characteristic dip of 7 seen in Figure 6(b) is still visible in Figure 7(c). 
Appendix B addresses the dependence of the empirical skewness on all method parameters. 

We also note that the estimates for the linear decay rate, Ku and K"acf are more accurate 
than the the estimates for the nonlinearity but by their nature they cannot distinguish the 
time series generated by the saddle-node normal form from a linear surrogate. 

The summative conclusion from Figure 7 is that one can detect the underlying nonlin- 
earity of the deterministic part of Equation (3.1) by observing Cemp, N2 or the skewness 7 
if the time series is moderately long. One can expect this to be true in the more general 
case of a time series generated by a system with deterministic dynamics close to a fold and 
additive noise. While one can in principle recover the underlying parameter n from the es- 
timates for Jfy, A^2 and a (see Thompson and Sieber [12]), the uncertainty in N2 propagates 
dramatically into uncertainty for /i (one has to divide by A'2)- 

The results of Figure 7 show the stationary case (e = 0). For slowly drifting time se- 
ries the same analysis will then have to be applied to parts of the time series in sliding 
windows. We demonstrate this for two geological times series in Section 5. For time series 
with rapidly drifting system parameter e the estimates for the softening, Cemp. N2 and the 
skewness 7, will not give a detectable difference from the corresponding linear time series 
before tipping. The conclusion one would draw from this absence of detectable softening 
would be correct: tipping happens when the hnear decay rate k reaches the critical value 
(or sUghtly later, see Thompson and Sieber [12]). The same holds if the noise level is low, 
which is equivalent to rapid drift (see Thompson and Sieber [12] for the transformation). 
In practice one uses the argument the other way round: if the nonlinear estimates do not 
give a value significantly different from zero one will conclude that the ratio between noise 
level and drifting speed is small. 
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5. Studies of Geological Time Series 

We now apply our nonlinear investigations to two ancient climate tippings. The first is the 
major transition marking the end of the latest ice age which occurred about 17,000 years 
before the present (yrs BP). The data used for this is a temperature reconstruction from 
the Vostok ice core deuterium record (Petit et al. [22]). The second more recent tipping is 
the ending of the Younger Dryas using the grayscale from the Cariaco basin sediments in 
Venezuela (Hughen [23]). This Younger Dryas event was a curious cooUng, about 1 1,500 
years ago, just as the Earth was warming up after the last ice age. It ended in a dramatic 
tipping point, when the Arctic warmed by 7°C in 50 years. This sudden ending has been 
related (Houghton [24]) to a switching-on of the global oceanic thermohaline circulation 
(THC). This switch-on is known from extensive theoretical studies (Dijkstra and Weijer 
[25]) to be at a saddle-node fold arising as a perturbation of a sub-critical pitchfork (see, 
for example, Rahmstorf [26] and Thompson and Sieber [12]). 




Figure 8: Two predicted potential energy suii'aces for (a) tlie end of tlie last glaeiation, using the Vostok ice- 
core record, and (b) the end of the Younger Dryas event when the Arctic warmed by 7°C in 50 years, using the 
grey-scale of basin sediment in Cariaco, Venezuela. The colour shows the deviation from a (time-dependent) 
parabola fitted to Uemp{x,t), three of which are illustrated above. Red signifies a positive deviation, blue a 
negative deviation. Time is given in years before the present (BP). The plotted potential was obtained via 
fcmp = ^fcmploglPemp) /2 and (3.7) and sliding windows of half the length of the record. 

Under linear time-series analysis of the directly preceding data, neither of these two 
events shows a strong trend in the stability propagator (see Figure 9). Meanwhile sample 
estimates of the underlying potential functions are shown in Figure 8. 

In Figure 8(a), for the end of the last glaeiation, we see clearly that the well is soften- 
ing (falling beneath the parabolic fit) on the high-temperature end while hardening (rising 
above the parabola) on the low-temperature end. So there is a strong nonlinear signal that 
a jump to higher temperatures may be pending (as indeed it was). We show that this signal 
is statistically significant in our full analysis in Figure 9. By comparison, the study of the 
Younger Dryas, illustrated in Figure 8(b), provided no significant nonlinear conclusions. 

For Figure 9 we detrended both time series using Gaussian filtering (see Figures 9(c) 
and (d)) and estimated the linear indicators Kacf and We observe weak to non-existent 
trends, leaving the evidence inconclusive at the linear level (grey and black curves Fig- 
ure 9(e) and (f)). The ratio between Kacf and Kjj, both shown in Figure 9(e) and (f), gives 
an estimate of the variance of the noise input. The unit for Kacf is units ofx per time 
step (so, for example for Fig. 9 (e) °C/Af), and, correspondingly, the unit for Ku is Af/°C 
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Figure 9: Linear and nonlinear coefficients for two paleo-climate time series. (a,c,e,g): End of last glaciation 
(snapshot of data from Petit et al. [22]); (b,d,f,h): End of Younger Dryas (snapshot of data from Hughen [23]). 
Only the black part of the data in panels (a) and (b) was used in the analysis. Panel (c) and (d) show the time series 
after detrending. Panel (e) and (f) show the linear indicators KXcF and Ku . Panel (g) and (h) show the means of the 
estimates Ccrap, ^2 and the skewness 7, compared to zero (thin vertical black line) and a histogram of estimates 



sampled from 500 random linear time series generated by the linear process (4.2) with 
window length w = N/2. 



2VM = icacf- Sliding 



(unit of Kacf divided by unit for o^). Both units are arbitrary such that we do not indicate 
the scale of the y-axis in Figure 9(e,f). 

[Lenton et al. this volume] discuss the linear analysis of both time series in greater 
detail. Their time series Vostok corresponds to Figure 9(a), their time series Cariaco cor- 
responds to Figure 9(b). (Note that [Lenton et al. this volume] use the deuterium proxy 
directly whereas Figure 9(a) shows the temperature reconstruction.) Specifically, Figures 2 
and 4, and Figures 9 and Al in [Lenton et al. this volume] discuss how the results depend 
on method parameters such as detrending bandwidth and sliding window. Their analysis 
finds that the presence of a significant trend depends strongly on both parameters. See 
[Lenton et al. this volume] , Table 1 (rows Vostok and Cariaco), for a summary of the 
evidence at the linear level. 

At the nonlinear level, the time series for the end of the last glaciation has a strong 
quadratic nonlinearity in dxU{x) where the potential well U is presumed to govern the 
deterministic part of the dynamics (Figure 9(g)). The indicators Cemp, N2 and the skew- 
ness 7 in Figure 9(g) are all far removed from what can be expected by chance in a linear 
time series. The histogram in the background of Figure 9(g,h) has been sampled from 500 
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random linear time series generated by the linear process given in Equation (4.2) with 
— 2\/ju = ktacf where JCacf is taken from the estimate shown in Figure 9(e). The percent- 
ages at the top of Figure 9(g,h) express how far in the tail of the histogram the mean of 
the quantity extracted from the geological time series is. We warn that the three quantities 
7, A'2 and Cgmp are not really three independent indicators as they all depend on the same 
estimate of the empirical density pcmp- 

Visual inspection of the well shape in Figure 8 confirms that the well is softening 
(bending downward) on the high-temperature end and hardening (bending upward) on the 
low-temperature end. So, at the nonlinear level the time series data close to the bottom of 
the well gives already evidence for a propensity to escape toward larger temperatures. The 
time series for the end of the Younger Dryas does not show evidence for strong nonlinearity 
of the underlying dynamics at the second-order level (Figure 9(h)). The values for Ccmp> N2 
and the skewness 7 can all be explained by randomness as the histograms of estimates for 
the 500 linear time series (as the histograms in Figure 9(h) show). Note that we did not 
include the previous transition (visible at the left end in Figure 9(b)) into our analysis as 
our aim was to infer the nonlinearity exclusively from the data near the tipping equilibrium. 

6. Conclusion 

The main message of the paper is that the linear analysis investigated by [Lenton et al. 
this volume] on its own, even if all estimates are accurate, does not contain all information 
necessary to estimate the probability of tipping over time. The result of the linear analy- 
sis is an estimate of the decay rate JC, which is taken relative to the time spacing of the 
measurements. Even if this decay rate has an identifiable trend to zero the probability of 
tipping over time is determined by the dominant nonlinear coefficient A'2 in conjunction 
with K and the noise level. 

We found that the dominant nonlinear coefficient A^2 is much harder to estimate ac- 
curately, so we also looked at proxies that indicate nonlinearity such as the skewness 7 
(as proposed by Guttal and Jayaprakash [19]) or the coefficient Cemp from Equation (3.5) 
(which is also a scalar signed measure for deviation from linear theory). Our study of the 
saddle-node normal form suggests that both proxies, 7 and Ccmp> give values that are easier 
to distinguish from random chance than the estimate for itself. A significantly non-zero 
value of either of these proxies indicates that a quadratic nonlinear term is present and 
gives its sign. 

However, we found the uncertainty for moderately long time series (A^ = 2000) too 
large to translate the proxy back into a quantitatively reliable estimate for the normal form 
parameter (this would be necessary to read off the probability for tipping from the tables 
in Thompson and Sieber [12]). We also found that the nonUnear proxies do not have dis- 
cernible trends when one approaches tipping, so it makes sense only to extract their overall 
mean from the time series but not the trend. 
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Appendix A. Sensitivity analysis for estimates of nonlinear quantities 

in palaeoclimate records 

Figure 10 shows how the results in Figure 9(g,h) depend on our method parameters. Specif- 
ically, we vary the length of the sliding window and the bandwidth of the Gaussian kernel 
when fitting the quasi-equiUbrium drift. In all panels the x-axis is the length of the slid- 
ing window relative to the overall length of the time series on a logarithmic scale (overall 
length is N = 525 for the end of the last glaciation, = 2646 for the end of the Younger 
Dryas). For example, the x-value of —1 corresponds to a window length w = N/2. The 
y-axis is the kernel bandwidth for the Gaussian fitting kernel relative to the overall length 
of the time series on a logarithmic scale (this enters the function kdeld as a resolution 
parameter). For example, the j'-value of —3 means that the resolution parameter of kdeld 
is 2^ = 8. 
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For each measured quantity we compare the mean of its value over all sliding win- 
dows to the percentiles of the linear surrogate. In this way contour levels close to or 100 
correspond to significant deviations from a time series generated by a linear process. The 
distributions of these time series are shown as gray background in Fig. 9(g,h). The value 
shown as a black vertical line in Figure 9(g,h) was obtained using the method parameters 
highlighted by a black circle in Figure 10. 

Fig. 10 gives evidence that the nonlinear features of the time series for the end of the 
last glaciation are robust as long as the detrending remains coarse grained. The end of the 
Younger Dryas does not show any nonlinear feature that is significantly different from the 
surrogates for any choise of method parameters. 



Appendix B. Dependence of skewness on distance to tipping point 

The observations in Figure 6(b) and Figure 7(c) appear to contradict the statement by Gut- 
tal and Jayaprakash [19] that skewness increases as one approaches tipping. Furthermore, 
there is a quantiative difference between Figure 6(b) and Figure 7(c), even though they 
show the same qualitative feature: non-monotonicity of the skewness in the control param- 
eter jU. Note that Guttal and Jayaprakash [19] studied the same type of tipping as Figure 6, 
however, not in the fold normal form but in a particular ecological model. In order to inves- 
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tigate this apparent contradiction we look at how the skewness depends on two parameters: 
the control parameter pL and the value at which we consider a realisation as escaped. More 
precisely, both figures, 6(b) and 7(c), study 



(Bl) 



where a = 1 without loss of generality and /i is constant. Any realisation of (B 1) escapes 
to — oo in finite time almost surely such that (B 1) does not possess a stationary density. In 
the calculation for Figure 6 we modified the process (B 1) by re-injecting the realisation at 
+00 whenever it escaped to —oa. In this way the process has a stationary density given by 
the Fokker-Planck equation (3.3) with —dJJ{x) = jj. —x^. However, this is not accurately 
describing the density that one is interested in: the probability density of the realisations 
under the condition that the realisation has not yet 'escaped'. Clearly, this conditional prob- 
ability depends on what one means by escaped. So, we introduce the escape boundary b as 
an additional parameter. Whenever, a realisation gets smaller than the value — y/JT — b (that 
is, it has run away beyond the local maximum of the potential by distance b) we stop the 
integration of (B 1). When one considers the set of all realisations before escape then this 
also leads to a stationary density. This was done for Figure 7 using b — — I + ^/JI (such 
that we count a realisation as escaped whenever it reaches —1). The re-injection used for 
Figure 6 and the choice of x = — 1 as escape indicator lead to the systematic difference 
in skewness 7 between the figures 7(c) and 6(b). In practice one encounters only the con- 
ditional probability as shown in Figure 7. So, it is important to check how the skewness 
depends on the somewhat arbitrary choice of the the value b indicating escape. Moreover, 
the estimate shown in Figure 7(c) suffers from the shorteness of the time series for jj. close 
to (which is typical in practice). 
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Figure 1 1 : Skewness of probability density of non-escaped realisations 



Figure 1 1 shows the result of a more precise calculation: we simulate (B 1) for a large 
ensemble of realisations (« = 100000). At every time step t/^ we collect the «esc.A: reali- 
sations that have reached a value x < —y/Ji—b. They count as escaped and have to be 
discarded. In order to avoid depletion of the ensemble we start n^sc,k new realisations at 
For each newly started realisation we pick the current value of one randomly chosen 
non-escaped realisation as the inital value. This process reaches a stationary density, which 
is an approximation of the conditional probability density of x under the condition that x 
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has not yet escaped beyond — b. We observe in Figure 1 1 that the depth of the dip 
in the skewness depends on b (indicated by grayscale) but that the location of the dip is 
uniformly larger than the critical value of the control parameter ji. 

It is possible that the simulations by Guttal and Jayaprakash [19] stayed consistently to 
the right of the curve of extreme skewness shown in Figure 1 1 such that only the increase 
is observed. Note that in their model the notion of skewness has the opposite sign, so one 
should consider the modulus of the skewness. 
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